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ABSTRACT 

We measure the 2-1 cumulant correlator power spectrum Cf 1 , a degenerate bispectrum, from the 
second data release of the Wilkinson Microwave Anisotropy Probe (WMAP). Our high resolution 
measurements with SpICE span a large configuration space (~ 168 x 999) corresponding to the possible 
cross-correlations of the maps recorded by the different differencing assemblies. We present a novel 
method to recover the eigenmodes of the correspondingly large Monte Carlo covariance matrix. We 
examine its eigenvalue spectrum and use random matrix theory to show that the off diagonal terms are 
dominated by noise. We minimize the % 2 to obtain constraints for the non-linear coupling parameter 
/jvl = 22 ±52(L7). 

Subject headings: cosmic microwave background — cosmology: theory — methods: statistical 
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1. INTRODUCTION 

Quantifying the non-Gaussianity in the cosmic mi- 
crowave background (CMB) puts constraints on infla- 
tionary models and possibly identifies non-linear effects 
from large scale structure. In addition, systematics and 
foregrounds might also produce non-Gaussian signatures; 
this possibly weakens constraints on the primordial and 
secondary non-Gaussianities. 

A natural phenomenological parametrization of non- 
Gaussian models is in terms of the perturbative non- 
linear coupling in the prim ordial curvature perturbation 
ijKomatsu fc Spergelll200l|) : 

d>(x) = $ L (x) + f NL ($|(x) - ($|(x)» , (1) 

where $l(x) denotes the linear Gaussian part of the 
Bardeen curvature and /jvl is the non- linear coupling 
parameter. The resulting leading-order non-Gaussianity 
is at the three-poi nt level. Thus, the thre e point correla- 
tion function (e.g.. lChen fc Szapudill2005l and references 
therein) or its spherical harmonic transform, the bispec- 
trum, directly estimate the leading order effect. 

The bispectrum has been used extens i vely for studying 
non-G aussianity (iKomatsu et al Jl2005t iC reminel li et al.l 
[20061 iMedeiros fc Contaldil 120061 iLiguori et alJ 120061: 
Cabe lla et all I2006L In previous measurements, the 
pseudo-bispectrum was used, which, like the pseudo-GYs, 
ignores in detail the effects of the complicated geometry 
induced by Galactic cut and cut-out holes. Pixel space 
statistics, such as the three-point correlation function, 
trivially deconvolve the geometric effects, as the convo- 
lution kernel is diago nal in pixel space. Indeed, SpICE 
(jSzapudi et al.ll2001al Spatially Inhomogenous Correla- 
tion Estimator) uses this simple fact to estimate the an- 
gular power spectrum withou t explicitly inverting the 
Mui kernel (Hivon et al.ll2002j) . For the bispectrum, the 
convolution kernel is even more complex than for the 
Cj's, therefore pixel space methods are advantagous. We 
use the fact that the SpICE algorithm can be used to 
calculate (deconvolv ed) power spectrum o f cumulant cor- 
relators (6T N 5T M ) llSzanudi et alJll99l . These are de- 
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generate N + M-point correlation functions, and their 
power spectra correspond to integrated N + M — 1 poly- 
spectra. In this paper we focus on the 2-1 cumulant 
correlator power spectrum which is directly related to 
bispectrum configurations. These contain less configra- 
tion information than the f ull bispectrum, althou gh re- 
tain more than the skewness l|Komatsu et al.l2~()">; . Note 
that in terms of Jnl, with near optimal weighting, the 
statistical power of t he skewness is nearly as optimal as 
the f ull bispectrum l)Komatsu et aLlfeOQSt iSpergel et alJ 
2006). For models with non-trivial configuration depen- 
dence, this is not the case. 

We introduce the power spectrum of the 2-1 cumu- 
lant correlator together with the corresponding theoret- 
ical theoretical predictions. Statistical estimation in the 
data and simulations, and the theoretical calculations are 
are described in §3. In §4 we investigate the covariance 
matrix of the measurements, the x 2 analysis of Jnl- We 
summarize and discuss the work in §5. 

2. POWER SPECTRUM OF 2-1 CUMULANT CORRELATOR 

For the temperature fluctuation field T(n) in CMB, 
the 2-1 cumulant correlator is simply expressed as 
(T(n) 2 T(m)), where the ensemble average can be re- 
placed by spatial (angular) average due to the assumed 
rotational invariance and ergodicity of the universe. Its 
Fourier or spherical harmonic transform corresponds to 
a set of summ ed (integrated) bispectrum configurations 
<ICooravll206"l : 



C? 1 




l)(2/ 2 + l) 



4tt(2/ + 1) 



(2) 

Here -Bz^z is the bispectrum and Wi is the multiple of 
the pixel window function and beam function of the CMB 
map. 

The above equation can be used to turn a theoreti- 
cal prediction for the full bispectrum into a prediction 
for the 2-1 cumulant correlator pow er spectrum. We fol- 
low cl osely the method described in IKomatsu fc Snergel 
l)2001|) to predict the full bispectrum using our modified 
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version of CAMB 2 . Then we sum the above Equation [3] 
in h and h up to k = 2000 where Wi ~ 0.001. Because 
of the linear dependency of the bispectrum on /nl, we 
perform the calculation with /jvi, = 1. 

Due to the similarity of cumulant correlators to two- 
point correlation functions, the analogous technique can 
be used for them as for measuring Cj. We use cross cor- 
relations of a triplet of maps for each measurement to 
avoid the uncertaintie s in the noise bias. The SpICE 
ijSzapudi et al.l l2001aj) algorithm uses harmonic trans- 
forms to calculate fast corr elation functions of tw o maps, 
and Legendre integration ijSzapudi et alJl2001b|) to ob- 
tain the final power spectrum. To obtain 2-1 cumulant 
correlators, the first two maps have to be multiplied first, 
as it is described next in detail. The final power spectrum 
can be directly compared with the theoretical prediction, 
as it is fully corrected for the complicated pixel geome- 
try of the underlying maps. Our measurement is the 
first such bispectrum measurement, where the geometry 
is accurately taken into account. 

3. MEASUREMENTS AND PREDICTIONS 

We use WMAP three y ear coadded foreground reduced 
sky maps i|Jarosik et al.ll2006L hereafter WMAP2) 3 with 
the resolution of n s ide = 512. There are 8 differencing as- 
semblies differencing assemblies (DAs) for the Q, V, and 
W bands. We calculated the cross correlations among 
these 8 maps and there are (8 x 7)/2 x 6 = 168 combi- 
nations. The factor of 1/2 is explained by the invariance 
of the 2-1 cumulant correlator under the exchange of the 
first two maps. We index the triplet of maps used in the 
cross-correlations by X and denote the power spectrum 
by Cf (X). First we prepare a "square temperature" 
map by multiplying pixel by pixel the first two maps 
in a triplet. The resulting map is cross-correlated by 
the third map with SpICE. We use the more conserva- 
tive KpO mask since non-Gaussianities are expected to be 
more sensitive to foregrounds than the power spectrum. 
It takes about 330 minutes to calculate all 168 spectra on 
a 2 GHz class CPU for the WMAP data (or for one set of 
simulations). Figure 1 displays typical measurements to- 
gether with theoretical predictions and simulation results 
for the triplet X=(Q1,Q2,V1). 

Qualitatively, there is no obvious sign of non- 
Gaussianities. To obtain accurate constraints, we es- 
timate covariances from a set of Monte Carlo simula- 
tions and fit the /jvl parameter. Gaussian simulations 
(/jvx=0) suffice, since the covariance i s dominated by 
the G aussian noise as long as /ml < 500 ( Ko matsu et al.l 
2003), which we already know to be true. We generated 
200 simulations with SYNFAST in HEALPix package 4 . 
We used the WMAP first year (WMAP1) power spec- 
trum available from the Lambda website 5 for input. It 
is the best fit ACDM model using a scale-dependent (run- 
ning) primordial spectral index, using WMAP1, CBI and 
ACBAR CMB data, plus the 2dF and Lyman-alpha data. 
For each simulat ion we generate 8 D A maps closely mim- 
icking the data ijSpergel et alJl2003|) . These maps repre- 
sent the same realization of the CMB but with different 

2 http://camb.info/ 

3 http://lambda.gsfc.nasa.gov/ 

4 http://healpix.jpl.nasa.gov/ 

5 The WMAP2 power spectrum was not available at this writing 
but the difference should be insignificant for our purposes. 
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Fig. 1. — A typical results of the 2-1 cumulant power spectrum 
normalized to the triplet X = (Ql, Q2, VI). The lower solid curve 
is the theoretical prediction for /nl = —100; the upper solid line 
is the standard deviation from 200 simulations; the dotted line 
corresponds to the absolute values from WMAP2. Inset: the data 
(without absolute value) and theory on a linear scale renormalized 
by the standard deviation measured from the simulations. 



noise and beam for each. Our noise model is somewhat 
simplistic, and uses Gaussian realizations with standard 
deviation <7o j\J N oos , where the effective number of ob- 
servations N b s varies across the sky and for different 
DAs, and ern is a constant for each DA. This does not 
take into account possible correlations in the noise, but 
at this writing no noise maps were available for WMAP2; 
our measurements in WMAP1 with more realistic noise 
indicate that the effects of the noise correlations are neg- 
ligible on our measurement (see the last section for de- 
tails). Finally, each simulation was analyzed exactly the 
same way as the data: we performed 33, 600 measure- 
ments in 16, 00 maps. 

4. THE COVARIANCE MATRIX AND x 2 ANALYSIS 

Let us extend our notation and label our measurements 
and predictions for triplet X as C 21 (A, s), where s stands 
for the theory (s = 0), one of the 200 simulations(s = 
1,...,200), or the WMAP data(s = 201). Our goal is 
to obtain quantitative constraints for j^L using these 
results. In particular, we focus on the covariance matrix 
(CM). 

4.1. The Covariance Matrix 

The simulations can be used to obtain an experimental 
CM the standard way. Without any further binning, the 
CM is a square matrix of size p — 168 x 999 = 167832 
(999 stands for I = 2, ...,1000). A matrix of this size 
can only be inverted with supercomputers. To overcome 
this problem, we introduce a novel, generally applicable 
technique which speeds up the inversion of experimental 
covariance matrices. We show that the singular value 
decomposition (SVD) of the CM is related to another 
matrix of smaller size governed by the number of simu- 
lations. 

Let M be the matrix where each raws and columns cor- 
respond to configurations and simulations, respectively. 
The SVD of this m x n matrix M (the number of rows 
m is typically larger than the number of columns n) 
is given by M = UAV T . Here U is a m x n column- 
orthonormal matrix, A is a n x n diagonal matrix with 
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non-negative elements, and V is a n x n orthonormal 
matrix. In this notation the CM can be expressed as 
l/nMM T . The corresponding SVD is U{l/nK 2 )U T . 
This can obtain more efficiently by solving the dual prob- 
lem M T M — V T A 2 V. Once the eigenvalues and V are 
obtained from this smaller matrix, U = MVA -1 can be 
calculated. This procedure is much faster than the direct 
calculation of the covariance matrix, as long asm> n. 
Also, it is clear from the arguments, that the CM has at 
most n non-degenerate eigenmodes. Unfortunately, us- 
ing only a small number of noisy eigenmodes means that 
we cannot fully exploit the information content of our 
data. While we have performed y 2 analysis using the 
full data set, a more reliable result can be obtained by 
compressing the data. 

4.2. Binning 

It would be desirable to combine different measure- 
ments of the same I with inverse variance weighting to 
compress our data set. Such a combi nation for the angu- 
lar po wer spectrum is nearly optimal ijFosalba fe SzapudU 
120041) . The generalization is a bit more complex, since 
the pixel and beam window functions for a particular 
triplet in Equation [21 cannot be decomposed, therefore 
one cannot simply obtain corrected Cf 1 estimates inde- 
pendent of window functions. 

To overcome this problem we introduce a new unbiased 
quantity, 



Cf"(X,s) 



(X,0) 



/nl — 1 



(3) 

the cumulant correlator power spectrum normalized ac- 
cording to the theory to a arbitrary triplet X — 
(Vl,V2,Ql) 6 at a fiducial value f NL = 1. In this 
normalization the window factors approximately cancel, 
therefore we can obtain an inverse variance weighted 



ET=iMX)C?HX,s) 



cf 1 



(s) 



(4) 



where w(x) is proportional to the inverse variance mea- 
sured in the simulations. The advantage of the in- 
verse variance weighting over the noisy covariance ma- 
trix is that each weight is determined with high accu- 
racy 1/V200 ~ 0.07, while, as we will see, modes of the 
matrix are significantly affected by the noise. Therefore, 
although our numerical method of the previous section 
allows us to handle a matrix of very large size, we opted 
for the inverse variance weighted estimators to suppress 
the noise more effectively. 

4.3. Random Matrix Theory 

To assess the level of noise in our final covariance ma- 
trix, we use some results from random matrix th eory. 
("e.g.. iSengupta fc Mitra1ll999t lLaloux et alJll999j) . As 
before, let us assume that a CM C is constructed as 
C = l/nMM T , from M, an to x n rectangular matrix. 
Let us further assume that M is composed of indepen- 
dent, identically distributed random variables with zero 
mean and unit variance. In the limit m — > oo, n — > oo 

6 Different choices of X make no difference for our results, so we 
will keep this normalization for the rest of the paper. 



while Q — n/m is kept fixed, the eigenvalue spectrum 
will tend to: 



Pc(X) 



Q 



A) (A — Xmin) 



2tt 



A 



\ m ax i 



1/Q±2y/T/Q, 



(5) 
(6) 



with A G [X m in, A„ 
Pc (A) is defined as 



„]. The eigenvalue density of C, 



PC (A) 



1 dn(X) 
to dX 



(7) 



where n(X) is the number of eigenvalues of C less than 
A. 

In our numerical construction of the covariance matrix, 
if C^were independent of each other, we would have a 
random matrix. How significant are the correlations? 
According to Figure 2, the density of eigenvalues of our 
CM appears to be consistent with the random matrix 
theory with Q = 999/200, with possibly a slight devia- 
tion for the largest eigenvalues. Comparing with random 
matrix simulations, we find that Q — 950/200 would pro- 
duce the same effect. While we keep Q — 999/200 in our 
analysis, this might be a sign of small correlations in the 
CM at the 5% level; this could effect our final % 2 and our 
final error-bars only slightly. 
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Fig. 2. — The eigenvalue spectrum for a set of power mapped 
CMs from 200 WMAP2 simulations and from the 200 random ma- 
trix simulations. The original matrix corresponds to q = 1. To 
display more detail, we plot unbinned (sorted) eigenvalues of our 
CM along with those obtained from random matrix simulations. 
The spectacular agreement indicates that the random matrix as- 
sumption is fairly accurate. 

To further investigate the randomness of the covari- 
ance matrix, we apply the technique of power mapping 
to our CM and random matrix simulations for compar- 
is on. The qth power m apping of a matrix C is defined 
as llGuhr fc Kalberil200l : 

Cl? = S ign(C kl )\C kl \ q - (8) 

iGuhr fc Kalberl l)2003fl show that the eigenvalue spectral 
density of power mapped correlation matrices can detect 
correlations otherwise buried in the noise. The variance 
of individual elements scales as 1/n 9 / 2 , i.e. effective num- 
ber of simulations increased to n q by power mapping. In 
our case, no significant deviation from random matrix 
theory appears as a result of power mapping (see Fig- 
ure 2) , which shows that the covariance is dominated by 
random noise. 
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5. RESULTS AND DISCUSSION 

Given that the CM of our measurement is consistent 
with the random matrix assumption at the 5% level, and 
that from 200 Monte Carlo simulations 7% fluctuations 
are expected for each of its elements, we conclude that 
it is consistent with our simulations to use a diagonal 
X 2 rather than weighting with the noisy off diagonal ele- 
ments of CM. We found that using diagonal x 2 is entirely 
robust when we vary our binning scheme or apply no bin- 
ning. On the contrary, using the noisy eigenmodes of the 
CM produces somewhat unstable results. 

Minimizing \ 2 as function of f^L gives /nl = 22 ± 
52 (ler), which is the final result of our paper. The min- 
imum reduced \ 2 1S about 1.1 for 998 degrees of free- 
dom. According to our discussion §4.3 there might be an 
additional 5% uncertainty on these results because we 
neglected off diagonal correlations. Our r esults are en- 
tirely consistent with ISpergel et al.l l|2006|) despite that 
we did not weight with the theory. The lack of opti- 
mal weighting is compensated by the increased config- 
uration dependence of our statistic. Note that we did 
not attempt to correct for point source contamination, 
but this should hav e a negligible effect on our results 
l|Spergel et al.ll200^ . 

To test the robustness of our constraints, we divided 
the 200 simulations to two sets of 100 simulations, then 
repeated the full statistical analysis with each set. We 



found consistent results: /jvl = 24 ± 54 and Jnl = 
21 ± 49, respectively. 

To test the degree to which noise correlations might af- 
fect our covariance, we also repeated o ur analysis of the 
CM w ith the WMAP1 simulations of iChen fc Szanudi 
(2005). where the correlated noise were used. The eigen- 
value spectrum is virtually identical to what we get for 
WMAP2 simulations. 

Finally, we repeated our analysis using the less conser- 
vative Kp2 mask and obtained Jml = —78 ± 52. This 
shows that the statistical variance of our results is al- 
ready of the same order of magnitude as the possible 
effect of foreground corrections near the galactic plane 
ijSpergel et al.ll2006|) . therefore it would be difficult to 
improve on significantly with the present data set. 
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